
# ----------------- Replication Code -----------------

# More Turnover, Less Turnout? Domestic Migration and Political Participation across Communities
# Giuliana Pardelli and Alexander Kustov (2024)
# British Journal of Political Science

# ----------------------------------------------------

library(lmtest)
library(sandwich)
library(multiwayvcov)
library(stargazer)
library(fixest)

library(maps)         
library(maptools)     
library(sp)          
library(RColorBrewer) 
library(classInt)

library(dplyr)
library(rgeos)
library(rgdal)

library(sf)
library(ggplot2)
library(ggthemes)
library(ggspatial)

library(panelr)
library(geepack)
library(modelsummary)

# ------------------------------------
#   Table 1 Cross-sectional Analysis 
# ------------------------------------

datamean <- read.csv("municipal-cross-section.csv")

ols1.avg <- lm(turnout_rate_fr.avg ~ scale(outmigrants.stock.pc) 
               + state
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor
               + latitude + longitude + distance_to_capital + distance_to_coast
               , data = datamean, weights = pop)

ols2.avg <- lm(turnout_rate_fr.avg ~ scale(inmigrants.stock.pc) 
               + state
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
               + latitude + longitude + distance_to_capital + distance_to_coast
               , data = datamean, weights = pop)

ols3.avg <- lm(turnout_rate_fr.avg ~ scale(turnover.stock.pc)
               + state
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor
               + latitude + longitude + distance_to_capital + distance_to_coast
               , data = datamean, weights = pop)

ols4.avg <- lm(turnout_rate_fr.avg ~ scale(net.mig.stock.pc)
               + state
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
               + latitude + longitude + distance_to_capital + distance_to_coast
               , data = datamean, weights = pop)

ols5.avg <- lm(turnout_rate_fr.avg ~ scale(outmigrants.stock.pc) + scale(inmigrants.stock.pc)
                + state
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                + latitude + longitude + distance_to_capital + distance_to_coast
                , data = datamean, weights = pop)


cl_ols1.avg <- coeftest(ols1.avg, vcov=function(x) cluster.vcov(x, datamean$state))[,2]  
cl_ols2.avg <- coeftest(ols2.avg, vcov=function(x) cluster.vcov(x, datamean$state))[,2] 
cl_ols3.avg <- coeftest(ols3.avg, vcov=function(x) cluster.vcov(x, datamean$state))[,2]
cl_ols4.avg <- coeftest(ols4.avg, vcov=function(x) cluster.vcov(x, datamean$state))[,2]
cl_ols5.avg <- coeftest(ols5.avg, vcov=function(x) cluster.vcov(x, datamean$state))[,2]

# Table 1 (left panel)
stargazer(ols1.avg, ols2.avg, ols5.avg, ols3.avg, ols4.avg,
          se=list(cl_ols1.avg, cl_ols2.avg, cl_ols5.avg, cl_ols3.avg, cl_ols4.avg),
          dep.var.caption = "",
          title="Cross-Sectional Analysis: Migration Stock Measure", 
          dep.var.labels=c("Average Turnout Rate (2000-2010)"),
          covariate.labels=c("Out-migration", "In-migration", "Migration Turnover", "Net Migration"),
          no.space=TRUE, column.sep.width = "10pt", 
          omit = c("Constant", "state", "electorate", "pop", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                   "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor",
                   "latitude", "longitude", "distance_to_coast", "distance_to_capital"),
          omit.stat=c("rsq", "ser", "f"), notes = "", notes.append = F,  notes.label = "",
          star.char = c("+", "*", "**", "***"), star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          add.lines = list(c("State FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE", "N/A", "N/A", "N/A", "N/A", "N/A"),
                           c("Control Variables", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Geographic controls", "Yes", "Yes", "Yes", "Yes", "Yes")))

# ----------------------------
#   Table 1 Panel Analysis
# ----------------------------

# Panel data with migration flow measures
mgp <- read.csv("municipal-yearly-panel.csv")
mgp = panel(mgp, ~ id_municipio + year)

m1t.r <- feols(turnout_rate_fr ~ scale(outmig.flow.pc) 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio + year, 
               data = mgp, vcov_cluster(~id_municipio), weights = ~pop)
summary(m1t.r)

m2t.r <- feols(turnout_rate_fr  ~ scale(inmig.flow.pc)  
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio + year, 
               data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m3t.r <- feols(turnout_rate_fr ~ scale(turnover.flow.pc) 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio + year, 
               data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m4t.r <- feols(turnout_rate_fr ~ scale(net.mig.flow.pc) 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio + year, 
               data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m5t.r <- feols(turnout_rate_fr ~ scale(outmig.flow.pc) + scale(inmig.flow.pc) 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio + year, 
               data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

# Table 1 (right panel)
etable(m1t.r, m2t.r, m5t.r, m3t.r, m4t.r,
       drop = c("pop", "local.election", "electorate", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor"),
       dict = c("turnout_rate_fr"="Biennial Turnout Rate (2002-2010)", "scale(outmig.flow.pc)" = "Out-migration", "scale(inmig.flow.pc)"="In-migration", "scale(turnover.flow.pc)"="Migration Turnover", "scale(net.mig.flow.pc)"="Net Migration"),
       title = "Panel Analysis: Migration Flow Measure",
       vcov = ~id_municipio, fitstat = ~ ar2 + n, 
       style.tex = style.tex("aer"), digits = "r3", tex = TRUE)


# ----------------------------------------
#       Appendix Tables and Figures 
# ----------------------------------------

## Maps ##

# Map: Municipalities
load("municipal-shapefile-2010.Rdata")
brazil2010_f <- fortify(mun)
brazil2010_f <- brazil2010_f %>%
  dplyr::rename(id_municipio = code_muni)
brazil2010_f <- left_join(brazil2010_f, datamean, by = "id_municipio")

# Map: Country
brazil <- st_read(dsn = "brazil_border/BRA_adm0.shp")
brazil <- st_transform(brazil, st_crs(brazil2010_f))

# --------------------------------
#     Figure A1: Turnout Rate
# --------------------------------

tiff("map_turnout.tiff", units="in", width=5, height=5, res=300)

# Average Turnout Across Brazilian Municipalities (2000-2010)
ggplot() +
  geom_sf(data = brazil2010_f, aes(fill = turnout_rate_fr.avg/100),color = NA) +
  scale_fill_gradient(low = "white", high = "#006699", 
                      name = "Turnout Rate") +
  theme_map() +
  ggtitle("Average Turnout Across Brazilian Municipalities (2000-2010)") +
  geom_sf(data = brazil, colour="gray70", fill = NA) +
  theme(legend.direction = "vertical",
        legend.justification = c("left", "bottom"),
        legend.spacing.y = unit(0, "mm"), 
        aspect.ratio = 1,
        legend.background = element_blank())
dev.off()


# -----------------
#    Figure A2 
# -----------------

# ------------------
#   In-Migration
# ------------------  

tiff("map_inmigration.tiff", units="in", width=5, height=5, res=300)

ggplot() +
  geom_sf(data = brazil2010_f, aes(fill = inmigrants.stock.pc), color = NA) +
  scale_fill_gradient(low = "white", high = "#006699", 
                      name = "In-migration Share") +
  theme_map() +
  #ggtitle("In-Migration") +
  geom_sf(data = brazil, colour="gray70", fill = NA) +
  theme(legend.direction = "vertical",
        legend.justification = c("left", "bottom"),
        legend.spacing.y = unit(0, "mm"), 
        aspect.ratio = 1,
        legend.background = element_blank())
dev.off()

# ------------------
#   Out-Migration
# ------------------  

tiff("map_outmigration.tiff", units="in", width=5, height=5, res=300)

ggplot() +
  geom_sf(data = brazil2010_f, aes(fill = outmigrants.stock.pc), color = NA) +
  scale_fill_gradient(low = "white", high = "#006699", 
                      name = "Out-migration Share") +
  theme_map() +
  #ggtitle("Out-Migration") +
  geom_sf(data = brazil, colour="gray70", fill = NA) +
  theme(legend.direction = "vertical",
        legend.justification = c("left", "bottom"),
        legend.spacing.y = unit(0, "mm"), 
        aspect.ratio = 1,
        legend.background = element_blank())
dev.off()

# ------------------
#   Net Migration
# ------------------ 

tiff("map_netmigration.tiff", units="in", width=5, height=5, res=300)

ggplot() +
  geom_sf(data = brazil2010_f, aes(fill = net.mig.stock.pc),color = NA) +
  scale_fill_gradient2(low = "red", mid = "white", high = "#006699", midpoint=0, 
                       name = 'Net Migration') +
  theme_map() +
  geom_sf(data = brazil, colour="gray70", fill = NA) +
  theme(legend.direction = "vertical",
        legend.justification = c("left", "bottom"),
        legend.spacing.y = unit(0, "mm"), 
        aspect.ratio = 1,
        legend.background = element_blank())
dev.off()

# ----------------------
#   Migration Turnover
# ---------------------- 

tiff("map_turnover.tiff", units="in", width=5, height=5, res=300)

ggplot() +
  geom_sf(data = brazil2010_f, aes(fill = turnover.stock.pc),color = NA) +
  scale_fill_gradient(low = "white", high = "#006699", 
                      name = "Migration Turnover") +
  theme_map() +
  #ggtitle("Migration Turnover") +
  geom_sf(data = brazil, colour="gray70", fill = NA) +
  theme(legend.direction = "vertical",
        legend.justification = c("left", "bottom"),
        legend.spacing.y = unit(0, "mm"), 
        aspect.ratio = 1,
        legend.background = element_blank())
dev.off()



# -----------------
#     Table A1 
# -----------------

ols1.avg.loc <- lm(turnout_rate_fr.local ~ scale(outmigrants.stock.pc) 
                   + state
                   + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                   + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor
                   + latitude + longitude + distance_to_capital + distance_to_coast
                   , data = datamean, weights = pop)

ols2.avg.loc <- lm(turnout_rate_fr.local ~ scale(inmigrants.stock.pc) 
                   + state
                   + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                   + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                   + latitude + longitude + distance_to_capital + distance_to_coast
                   , data = datamean, weights = pop)

ols3.avg.loc <- lm(turnout_rate_fr.local ~ scale(turnover.stock.pc)
                   + state
                   + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                   + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor
                   + latitude + longitude + distance_to_capital + distance_to_coast
                   , data = datamean, weights = pop)

ols4.avg.loc <- lm(turnout_rate_fr.local ~ scale(net.mig.stock.pc)
                   + state
                   + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                   + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                   + latitude + longitude + distance_to_capital + distance_to_coast
                   , data = datamean, weights = pop)

ols5.avg.loc <- lm(turnout_rate_fr.local ~ scale(outmigrants.stock.pc) + scale(inmigrants.stock.pc)
                    + state
                    + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                    + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                    + latitude + longitude + distance_to_capital + distance_to_coast
                    , data = datamean, weights = pop)


cl_ols1.avg.loc <- coeftest(ols1.avg.loc, vcov=function(x) cluster.vcov(x, datamean$state))[,2]  
cl_ols2.avg.loc <- coeftest(ols2.avg.loc, vcov=function(x) cluster.vcov(x, datamean$state))[,2] 
cl_ols3.avg.loc <- coeftest(ols3.avg.loc, vcov=function(x) cluster.vcov(x, datamean$state))[,2] 
cl_ols4.avg.loc <- coeftest(ols4.avg.loc, vcov=function(x) cluster.vcov(x, datamean$state))[,2] 
cl_ols5.avg.loc <- coeftest(ols5.avg.loc, vcov=function(x) cluster.vcov(x, datamean$state))[,2] 

# Table A1
stargazer(ols1.avg.loc, ols2.avg.loc, ols5.avg.loc, ols3.avg.loc, ols4.avg.loc,
          se=list(cl_ols1.avg.loc, cl_ols2.avg.loc, cl_ols5.avg.loc, cl_ols3.avg.loc, cl_ols4.avg.loc),
          dep.var.caption = "",
          title="Cross-Sectional Analysis: Average Migration Shares and Turnout in Local Elections (2000-2008)", 
          dep.var.labels=c("Average Turnout Rate in Local Elections"),
          covariate.labels=c("Out-migration", "In-migration", "Migration Turnover", "Net Migration"),
          no.space=TRUE, column.sep.width = "10pt", 
          omit = c("Constant", "state", "electorate", "pop", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                   "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor",
                   "latitude", "longitude", "distance_to_coast", "distance_to_capital"),
          omit.stat=c("rsq", "ser", "f"), notes = "", notes.append = F,  notes.label = "",
          star.char = c("+", "*", "**", "***"), star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          add.lines = list(c("State FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional Covariates", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Geographic controls", "Yes", "Yes", "Yes", "Yes", "Yes")))

# ------------
#   Table A2 
# ------------

ols1.avg.nat <- lm(turnout_rate_fr.national ~ scale(outmigrants.stock.pc) 
                   + state
                   + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                   + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor
                   + latitude + longitude + distance_to_capital + distance_to_coast
                   , data = datamean, weights = pop)

ols2.avg.nat <- lm(turnout_rate_fr.national ~ scale(inmigrants.stock.pc) 
                   + state
                   + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                   + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                   + latitude + longitude + distance_to_capital + distance_to_coast
                   , data = datamean, weights = pop)

ols3.avg.nat <- lm(turnout_rate_fr.national ~ scale(turnover.stock.pc)
                   + state
                   + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                   + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor
                   + latitude + longitude + distance_to_capital + distance_to_coast
                   , data = datamean, weights = pop)

ols4.avg.nat <- lm(turnout_rate_fr.national ~ scale(net.mig.stock.pc)
                   + state
                   + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                   + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                   + latitude + longitude + distance_to_capital + distance_to_coast
                   , data = datamean, weights = pop)

ols5.avg.nat <- lm(turnout_rate_fr.national ~ scale(outmigrants.stock.pc) + scale(inmigrants.stock.pc)
                    + state
                    + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                    + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                    + latitude + longitude + distance_to_capital + distance_to_coast
                    , data = datamean, weights = pop)


cl_ols1.avg.nat <- coeftest(ols1.avg.nat, vcov=function(x) cluster.vcov(x, datamean$state))[,2]  
cl_ols2.avg.nat <- coeftest(ols2.avg.nat, vcov=function(x) cluster.vcov(x, datamean$state))[,2] 
cl_ols3.avg.nat <- coeftest(ols3.avg.nat, vcov=function(x) cluster.vcov(x, datamean$state))[,2] 
cl_ols4.avg.nat <- coeftest(ols4.avg.nat, vcov=function(x) cluster.vcov(x, datamean$state))[,2] 
cl_ols5.avg.nat <- coeftest(ols5.avg.nat, vcov=function(x) cluster.vcov(x, datamean$state))[,2] 

#  Avg Turnout in National elections
stargazer(ols1.avg.nat, ols2.avg.nat, ols5.avg.nat, ols3.avg.nat, ols4.avg.nat,
          se=list(cl_ols1.avg.nat, cl_ols2.avg.nat, cl_ols5.avg.nat, cl_ols3.avg.nat, cl_ols4.avg.nat),
          dep.var.caption = "",
          title="Cross-Sectional Analysis: Average Migration Shares and Turnout in National Elections (2002-2010)", 
          dep.var.labels=c("Average Turnout Rate in National Elections"),
          covariate.labels=c("Out-migration", "In-migration", "Migration Turnover", "Net Migration"),
          no.space=TRUE, column.sep.width = "10pt", 
          omit = c("Constant", "state", "electorate", "pop", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                   "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor",
                   "latitude", "longitude", "distance_to_coast", "distance_to_capital"),
          omit.stat=c("rsq", "ser", "f"), notes = "", notes.append = F,  notes.label = "",
          star.char = c("+", "*", "**", "***"), star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          add.lines = list(c("State FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional Covariates", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Geographic controls", "Yes", "Yes", "Yes", "Yes", "Yes")))


# -------------
#   Table A3 
# -------------

# Panel A

olspc1.00 <- lm(turnout_rate_fr ~ scale(outmig.stock.pc) 
                + state
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor
                , data = mgp[which(mgp$year == 2000),], weights = pop)

olspc2.00 <- lm(turnout_rate_fr ~ scale(inmig.stock.pc) 
                + state
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                , data = mgp[which(mgp$year == 2000),], weights = pop)

olspc3.00 <- lm(turnout_rate_fr ~ scale(turnover.stock.pc) 
                + state
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                , data = mgp[which(mgp$year == 2000),], weights = pop)

olspc4.00 <- lm(turnout_rate_fr ~ scale(net.mig.stock.pc)
                + state
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                , data = mgp[which(mgp$year == 2000),], weights = pop)

olspc5.00 <- lm(turnout_rate_fr ~ scale(outmig.stock.pc) + scale(inmig.stock.pc) 
                + state
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                , data = mgp[which(mgp$year == 2000),], weights = pop)

cl_olspc1.00 <- coeftest(olspc1.00, vcov=function(x) cluster.vcov(x, mgp$state[which(mgp$year == 2000)]))[,2] 
cl_olspc2.00 <- coeftest(olspc2.00, vcov=function(x) cluster.vcov(x, mgp$state[which(mgp$year == 2000)]))[,2] 
cl_olspc3.00 <- coeftest(olspc3.00, vcov=function(x) cluster.vcov(x, mgp$state[which(mgp$year == 2000)]))[,2] 
cl_olspc4.00 <- coeftest(olspc4.00, vcov=function(x) cluster.vcov(x, mgp$state[which(mgp$year == 2000)]))[,2] 
cl_olspc5.00 <- coeftest(olspc5.00, vcov=function(x) cluster.vcov(x, mgp$state[which(mgp$year == 2000)]))[,2] 

# Panel B

olspc1.10 <- lm(turnout_rate_fr ~ scale(outmig.stock.pc) 
                + state
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor
                , data = mgp[which(mgp$year == 2010),], weights = pop)

olspc2.10 <- lm(turnout_rate_fr ~ scale(inmig.stock.pc) 
                + state
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                , data = mgp[which(mgp$year == 2010),], weights = pop)

olspc3.10 <- lm(turnout_rate_fr ~ scale(turnover.stock.pc) 
                + state
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                , data = mgp[which(mgp$year == 2010),], weights = pop)

olspc4.10 <- lm(turnout_rate_fr ~ scale(net.mig.stock.pc)
                + state
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                , data = mgp[which(mgp$year == 2010),], weights = pop)

olspc5.10 <- lm(turnout_rate_fr ~ scale(outmig.stock.pc) + scale(inmig.stock.pc) 
                + state
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                , data = mgp[which(mgp$year == 2010),], weights = pop)
summary(olspc5.10)

cl_olspc1.10 <- coeftest(olspc1.10, vcov=function(x) cluster.vcov(x, mgp$state[which(mgp$year == 2010)]))[,2] 
cl_olspc2.10 <- coeftest(olspc2.10, vcov=function(x) cluster.vcov(x, mgp$state[which(mgp$year == 2010)]))[,2] 
cl_olspc3.10 <- coeftest(olspc3.10, vcov=function(x) cluster.vcov(x, mgp$state[which(mgp$year == 2010)]))[,2] 
cl_olspc4.10 <- coeftest(olspc4.10, vcov=function(x) cluster.vcov(x, mgp$state[which(mgp$year == 2010)]))[,2] 
cl_olspc5.10 <- coeftest(olspc5.10, vcov=function(x) cluster.vcov(x, mgp$state[which(mgp$year == 2010)]))[,2] 

# Table A3 Panel A
stargazer(olspc1.00, olspc2.00, olspc5.00, olspc3.00, olspc4.00,
          se=list(cl_olspc1.00, cl_olspc2.00, cl_olspc5.00, cl_olspc3.00, cl_olspc4.00),
          dep.var.caption = "",
          title="Cross-Section: Migration Stock and Turnout (2000)", 
          dep.var.labels=c("Turnout Rate 2000, Local Elections"),
          covariate.labels=c("Out-migration", "In-migration", "Migration Turnover", "Net Migration"),
          no.space=TRUE, column.sep.width = "10pt", 
          omit = c("Constant", "state", "pop", "electorate", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                   "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor"),
          omit.stat=c("rsq", "ser", "f"), notes = "", notes.append = F,  notes.label = "",
          star.char = c("+", "*", "**", "***"), star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          add.lines = list(c("State FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional Covariates", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Geographic controls", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")))

# Table A3 Panel B
stargazer(olspc1.10, olspc2.10, olspc5.10, olspc3.10, olspc4.10,
          se=list(cl_olspc1.10, cl_olspc2.10, cl_olspc5.10, cl_olspc3.10, cl_olspc4.10),
          dep.var.caption = "",
          title="Cross-Section: Migration Stock and Turnout (2010)", 
          dep.var.labels=c("Turnout Rate 2010, National Elections"),
          covariate.labels=c("Out-migration", "In-migration", "Migration Turnover", "Net Migration"),
          no.space=TRUE, column.sep.width = "10pt", 
          omit = c("Constant", "state", "pop", "electorate", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                   "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor"),
          omit.stat=c("rsq", "ser", "f"), notes = "", notes.append = F,  notes.label = "",
          star.char = c("+", "*", "**", "***"), star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          add.lines = list(c("State FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional Covariates", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Geographic controls", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")))


# --------------------
#      Table A4
# --------------------

# Panel data with migration stock measures (1991, 2000, 2010) 
wv <- read.csv("municipal-census-waves-panel.csv")  

m1s.r <- feols(turnout_rate_fr ~ scale(outmigrants.pc)
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio + year, 
               data = wv, vcov_cluster(~id_municipio), weights = ~pop)

m2s.r <- feols(turnout_rate_fr  ~ scale(inmigrants.pc)
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio + year, 
               data = wv, vcov_cluster(~id_municipio), weights = ~pop)

m3s.r <- feols(turnout_rate_fr ~ scale(turnover.pc) 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio + year, 
               data = wv, vcov_cluster(~id_municipio), weights = ~pop)

m4s.r <- feols(turnout_rate_fr ~ scale(net.migration.pc) 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio + year, 
               data = wv, vcov_cluster(~id_municipio), weights = ~pop)

m5s.r <- feols(turnout_rate_fr ~ scale(outmigrants.pc) + scale(inmigrants.pc)
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio + year, 
               data = wv, vcov_cluster(~id_municipio), weights = ~pop)

# Table A4
etable(m1s.r, m2s.r, m5s.r, m3s.r, m4s.r,
       drop = c("pop", "electorate", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor"),
       dict = c("turnout_rate_fr"="Turnout Rate", "scale(outmig.flow.pc)" = "Out-migration", 
                "scale(inmig.flow.pc)"="In-migration", "scale(turnover.flow.pc)"="Migration Turnover", 
                "scale(net.mig.flow.pc)"="Net Migration",
                "scale(outmigrants.pc)" = "Out-migration", "scale(inmigrants.pc)" = "In-migration",
                "scale(turnover.pc)" = "Migration Turnover", "scale(net.migration.pc)" = "Net Migration"),
       vcov = ~id_municipio, fitstat = ~ ar2 + n, 
       title = "Two-Way Fixed Effects Analysis: Turnout Rate and Migration Shares (Stock Measure)",
       style.tex = style.tex("aer"), digits = "r3", tex = TRUE)


# ---------------
#    Table A5 
# ---------------

m1.int <- feols(turnout_rate_fr ~ scale(outmig.flow.pc)*local.election
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                | id_municipio, 
                data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m2.int <- feols(turnout_rate_fr ~ scale(inmig.flow.pc)*local.election  
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                | id_municipio, 
                data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m3.int <- feols(turnout_rate_fr ~ scale(turnover.flow.pc)*local.election
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                | id_municipio, 
                data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m4.int <- feols(turnout_rate_fr ~ scale(net.mig.flow.pc)*local.election 
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                | id_municipio, 
                data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m5.int <- feols(turnout_rate_fr ~ scale(outmig.flow.pc)*local.election 
                + scale(inmig.flow.pc)*local.election 
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                | id_municipio, 
                data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

# Table A5 
etable(m1.int, m2.int, m5.int, m3.int, m4.int,
       drop = c("pop", "electorate", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor"),
       dict = c("local.election" = "Local Elections", "turnout_rate_fr"="Turnout Rate", "scale(outmig.flow.pc)" = "Out-migration", "scale(inmig.flow.pc)"="In-migration", "scale(turnover.flow.pc)"="Migration Turnover", "scale(net.mig.flow.pc)"="Net Migration"),
       title = "Two-Way Fixed Effects Analysis: Turnout Rate and Migration Shares, Interaction with Local-Election Indicator",
       vcov = ~id_municipio, fitstat = ~ ar2 + n, 
       style.tex = style.tex("aer"), digits = "r3", tex = TRUE)

# ----------------
#    Table A6
# ----------------

library(spdep)
library(spgwr)
library(sphet)
library(geobr)
library(sf)
library(spatialreg)
library(lwgeom)
set.seed(1234)

# 2010 Census Data
map <- read_municipality(year = 2010)
map <- map %>%
  dplyr::rename(id_municipio = code_muni)

map <- left_join(map, mgp[which(mgp$year == 2010),], by = "id_municipio")

map$inmigrants.pc <- map$inmigrants.stock/map$pop
map$outmigrants.pc <- map$outmigrants.stock/map$pop
map$net.migration.stock <- map$inmigrants.stock - map$outmigrants.stock 
map$turnover.pc <- (map$outmigrants.stock + map$inmigrants.stock)/map$pop
map$net.migration.pc <- map$net.migration.stock/map$pop

map <- map[which(!is.na(map$pop) & !is.na(map$turnout_rate_fr) & !is.na(map$electorate) & !is.na(map$inmigrants.pc)),]

# Make geometries valid
map <- st_make_valid(map)
map <- st_buffer(map, dist = 0)

# spatial weights matrix
weights <- nb2listw(poly2nb(map, queen=T), style = "W", zero.policy=TRUE) 

# SAR Models
mod.sar1 <- lagsarlm(turnout_rate_fr ~ scale(inmigrants.pc) + scale(outmigrants.pc) 
                     + state 
                     + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                     + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                     , data = map, listw = weights, zero.policy=T, tol.solve=1e-12)

mod.sar2 <- lagsarlm(turnout_rate_fr ~ scale(turnover.pc)
                     + state 
                     + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                     + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                     , data = map, listw = weights, zero.policy=T, tol.solve=1e-12)

mod.sar3 <- lagsarlm(turnout_rate_fr ~ scale(net.migration.pc) 
                     + state 
                     + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                     + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                     , data = map, listw = weights, zero.policy=T, tol.solve=1e-12)

# 2000 Census Data
map00 <- read_municipality(year = 2000)

# Make geometries valid
map00$geom <- st_make_valid(map00$geom)

map00 <- map00 %>%
  dplyr::rename(id_municipio = code_muni)

# Merge with census data
map00 <- left_join(map00, mgp[which(mgp$year == 2000),], by = "id_municipio")

map00$inmigrants.pc <- map00$inmigrants.stock/map00$pop
map00$outmigrants.pc <- map00$outmigrants.stock/map00$pop
map00$net.migration.stock <- map00$inmigrants.stock - map00$outmigrants.stock 
map00$turnover.pc <- (map00$outmigrants.stock + map00$inmigrants.stock)/map00$pop
map00$net.migration.pc <- map00$net.migration.stock/map00$pop

map00 <- map00[which(!is.na(map00$pop) & !is.na(map00$turnout_rate_fr) & !is.na(map00$electorate)),]

# spatial weights matrix
weight00 <- nb2listw(poly2nb(map00, queen=T), style = "W", zero.policy=TRUE) 

# Models
mod.sar4 <- lagsarlm(turnout_rate_fr ~ scale(inmigrants.pc) + scale(outmigrants.pc) 
                     + state 
                     + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                     + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                     , data = map00, listw = weight00, zero.policy=T, tol.solve=1e-12)

mod.sar5 <- lagsarlm(turnout_rate_fr ~ scale(turnover.pc)
                     + state 
                     + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                     + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                     , data = map00, listw = weight00, zero.policy=T, tol.solve=1e-12)

mod.sar6 <- lagsarlm(turnout_rate_fr ~ scale(net.migration.pc) 
                     + state 
                     + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                     + log(pop_urb) + log(pop.college+1) + avg.capita.income + Gini + pc.poor 
                     , data = map00, listw = weight00, zero.policy=T, tol.solve=1e-12)

# Table A6 
stargazer(mod.sar4, mod.sar5, mod.sar6, mod.sar1, mod.sar2, mod.sar3,
          dep.var.caption = "",
          title="Spatial Autoregressive Analysis: Turnout Rate and Migration Shares (Stock Measure)", 
          dep.var.labels=c("Turnout Rate 2000, Local Election", "Turnout Rate 2010, National Election"),
          covariate.labels=c("Out-migration", "In-migration", "Migration Turnover", "Net Migration"),
          no.space=TRUE, column.sep.width = "10pt", 
          omit = c("Constant", "state", "pop", "electorate", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                   "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor"),
          keep.stat=c("rsq", "n"), notes = "", notes.append = F,  notes.label = "",
          star.char = c("+", "*", "**", "***"), star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          add.lines = list(c("State FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional Covariates", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")))

# -----------------
#     Table A7 
# -----------------

# ---------- Flow Measure ------------

rv.flow <- feols(turnout_rate_fr ~ scale(outmig.flow.pc) + scale(inmig.flow.pc) 
                 + log(electorate)
                 + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                 + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                 | id_municipio + year, 
                 data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

rv.flow2 <- feols(turnout_rate_fr ~ scale(turnover.flow.pc) 
                  + log(electorate)
                  + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                  + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                  | id_municipio + year, 
                  data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

vep.flow <- feols(turnout_rate_vep ~ scale(outmig.flow.pc) + scale(inmig.flow.pc)
                  + log(electorate)
                  + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                  + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                  | id_municipio + year, 
                  data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

vep.flow2 <- feols(turnout_rate_vep ~ scale(turnover.flow.pc) 
                   + log(electorate)
                   + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                   + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                   | id_municipio + year, 
                   data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

# ---------- Stock Measure ------------

rv.stock <- feols(turnout_rate_fr ~ scale(outmigrants.pc) + scale(inmigrants.pc)
                  + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                  + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                  | id_municipio + year, 
                  data = wv, vcov_cluster(~id_municipio), weights = ~pop)

rv.stock2 <- feols(turnout_rate_fr ~ scale(turnover.pc) 
                   + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                   + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                   | id_municipio + year, 
                   data = wv, vcov_cluster(~id_municipio), weights = ~pop)

vep.stock <- feols(turnout_rate_vep ~ scale(outmigrants.pc) + scale(inmigrants.pc)
                   + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                   + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                   | id_municipio + year, 
                   data = wv, vcov_cluster(~id_municipio), weights = ~pop)

vep.stock2 <- feols(turnout_rate_vep ~ scale(turnover.pc) 
                    + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                    + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                    | id_municipio + year, 
                    data = wv, vcov_cluster(~id_municipio), weights = ~pop)

# Table A7 
etable(rv.flow, rv.flow2, vep.flow, vep.flow2, rv.stock, rv.stock2, vep.stock, vep.stock2, 
       drop = c("pop", "electorate", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor"),
       dict = c("turnout_rate_fr"="Turnout Rate (RV)", "turnout_rate_vep" = "Turnout Rate (VEP)", 
                "scale(outmig.flow.pc)" = "Out-migration", 
                "scale(inmig.flow.pc)"="In-migration", 
                "scale(turnover.flow.pc)"="Migration Turnover", 
                "scale(outmigrants.pc)" = "Out-migration", 
                "scale(inmigrants.pc)"="In-migration", 
                "scale(turnover.pc)" = "Migration Turnover" ),
       title = "Panel Data Analysis: Turnout Rates (RV and VEP) and Migration Shares (Flow and Stock)",
       vcov = ~id_municipio, fitstat = ~ ar2 + n, 
       style.tex = style.tex("aer"), digits = "r3", tex = TRUE)

# ---------------
#    Table A8 
# ---------------

# ---------------- Flow Measures ----------------

m1t <- feols(log(turnout_firstround) ~ log(outmig_flow)
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
             + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
             | id_municipio + year, 
             data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m2t <- feols(log(turnout_firstround) ~ log(inmig_flow) + 
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
             + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
             | id_municipio + year, 
             data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m3t <- feols(log(turnout_firstround) ~ log(turnover.flow) 
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
             + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
             | id_municipio + year, 
             data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m4t <- feols(log(turnout_firstround) ~ scale(net.mig.flow) # Using the log excludes observations ≤ 0
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
             + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
             | id_municipio + year, 
             data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m5t <- feols(log(turnout_firstround) ~ log(outmig_flow) + log(inmig_flow)
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
             + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
             | id_municipio + year, 
             data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

# ------------ Stock Measures -------------

m1s <- feols(log(turnout_firstround) ~ log(outmigrants.stock)
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
             + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
             | id_municipio + year, 
             data = wv, vcov_cluster(~id_municipio), weights = ~pop)

m2s <- feols(log(turnout_firstround)  ~ log(inmigrants.stock)
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
             + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
             | id_municipio + year, 
             data = wv, vcov_cluster(~id_municipio), weights = ~pop)

m3s <- feols(log(turnout_firstround) ~ log(turnover.stock) 
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
             + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
             | id_municipio + year, 
             data = wv, vcov_cluster(~id_municipio), weights = ~pop)

m4s <- feols(log(turnout_firstround) ~ scale(net.migration.stock) # Using the log excludes observations ≤ 0
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
             + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
             | id_municipio + year, 
             data = wv, vcov_cluster(~id_municipio), weights = ~pop)

m5s <- feols(log(turnout_firstround) ~ log(outmigrants.stock) + log(inmigrants.stock)
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
             + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
             | id_municipio + year, 
             data = wv, vcov_cluster(~id_municipio), weights = ~pop)


# Table A8 
etable(m5t, m3t, m4t, m5s, m3s, m4s,
       drop = c("pop", "electorate", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor"),
       vcov = ~id_municipio, fitstat = ~ ar2 + n, 
       dict = c("log(turnout_firstround)"="Turnout Count (log)", "log(outmig_flow)" = "Out-migration (log)", 
                "log(inmig_flow)"="In-migration (log)", "log(turnover.flow)"="Migration Turnover (log)", 
                "scale(net.mig.flow)"="Net Migration",
                "log(outmigrants.stock)" = "Out-migration (log)", "log(inmigrants.stock)" = "In-migration (log)",
                "log(turnover.stock)" = "Migration Turnover (log)", "scale(net.migration.stock)" = "Net Migration"),
       title = "Panel Data Analysis: Turnout and Migration Counts (Flow and Stock)",
       style.tex = style.tex("aer"), digits = "r3", tex = TRUE)


# ------------
#  Table A9 
# ------------

# MCA-level: Panel Data (2000, 2010)
wva <- read.csv("mca-panel.csv")

wva$outmigrants.pc <- wva$outmigrants / wva$pop.tot
wva$inmigrants.pc <- wva$inmigrants / wva$pop.tot
wva$turnover.pc <- (wva$inmigrants + wva$outmigrants) / wva$pop.tot
wva$net.migration.pc <- (wva$inmigrants - wva$outmigrants) / wva$pop.tot

m1s.amc <- feols(turnout_rate_fr ~ scale(outmigrants.pc)
                 + log(electorate) + log(pop.tot) + log(people_between_18_and_70) + log(people_above_16) + log(area_km2)
                 + log(pop.urb) + log(pop_higher_educ) + avg.capita.income + Gini + pc.poor 
                 | id_amc + year, 
                 data = wva, vcov_cluster(~id_amc), weights = ~pop.tot)

m2s.amc <- feols(turnout_rate_fr  ~ scale(inmigrants.pc)
                 + log(electorate) + log(pop.tot) + log(people_between_18_and_70) + log(people_above_16) + log(area_km2)
                 + log(pop.urb) + log(pop_higher_educ) + avg.capita.income + Gini + pc.poor 
                 | id_amc + year, 
                 data = wva, vcov_cluster(~id_amc), weights = ~pop.tot)

m3s.amc <- feols(turnout_rate_fr ~ scale(turnover.pc) 
                 + log(electorate) + log(pop.tot) + log(people_between_18_and_70) + log(people_above_16) + log(area_km2)
                 + log(pop.urb) + log(pop_higher_educ) + avg.capita.income + Gini + pc.poor 
                 | id_amc + year, 
                 data = wva, vcov_cluster(~id_amc), weights = ~pop.tot)

m4s.amc <- feols(turnout_rate_fr ~ scale(net.migration.pc) 
                 + log(electorate) + log(pop.tot) + log(people_between_18_and_70) + log(people_above_16) + log(area_km2)
                 + log(pop.urb) + log(pop_higher_educ) + avg.capita.income + Gini + pc.poor 
                 | id_amc + year, 
                 data = wva, vcov_cluster(~id_amc), weights = ~pop.tot)

m5s.amc <- feols(turnout_rate_fr ~ scale(outmigrants.pc) + scale(inmigrants.pc)
                 + log(electorate) + log(pop.tot) + log(people_between_18_and_70) + log(people_above_16) + log(area_km2)
                 + log(pop.urb) + log(pop_higher_educ) + avg.capita.income + Gini + pc.poor 
                 | id_amc + year, 
                 data = wva, vcov_cluster(~id_amc), weights = ~pop.tot)

# Micro-Region level: Panel Data (2000, 2010)
wvm <- read.csv("microregion-panel.csv")

wvm$inmigrants.pc <- wvm$inmigrants/wvm$pop.tot
wvm$outmigrants.pc <- wvm$outmigrants/wvm$pop.tot
wvm$turnover.pc <- (wvm$outmigrants + wvm$inmigrants)/wvm$pop.tot
wvm$net.migration.pc <- wvm$net.migration/wvm$pop.tot

m1s.mcr <- feols(turnout_rate_fr ~ scale(outmigrants.pc)
                 + log(electorate) + log(pop.tot) + log(people_between_18_and_70) + log(people_above_16) + log(area_km2)
                 + log(pop.urb) + log(pop_higher_educ) + avg.capita.income + Gini + pc.poor 
                 | microregion + year, 
                 data = wvm, vcov_cluster(~microregion), weights = ~pop.tot)

m2s.mcr <- feols(turnout_rate_fr  ~ scale(inmigrants.pc)
                 + log(electorate) + log(pop.tot) + log(people_between_18_and_70) + log(people_above_16) + log(area_km2)
                 + log(pop.urb) + log(pop_higher_educ) + avg.capita.income + Gini + pc.poor 
                 | microregion + year, 
                 data = wvm, vcov_cluster(~microregion), weights = ~pop.tot)

m3s.mcr <- feols(turnout_rate_fr ~ scale(turnover.pc) 
                 + log(electorate) + log(pop.tot) + log(people_between_18_and_70) + log(people_above_16) + log(area_km2)
                 + log(pop.urb) + log(pop_higher_educ) + avg.capita.income + Gini + pc.poor 
                 | microregion + year, 
                 data = wvm, vcov_cluster(~microregion), weights = ~pop.tot)

m4s.mcr <- feols(turnout_rate_fr ~ scale(net.migration.pc) 
                 + log(electorate) + log(pop.tot) + log(people_between_18_and_70) + log(people_above_16) + log(area_km2)
                 + log(pop.urb) + log(pop_higher_educ) + avg.capita.income + Gini + pc.poor 
                 | microregion + year, 
                 data = wvm, vcov_cluster(~microregion), weights = ~pop.tot)

m5s.mcr <- feols(turnout_rate_fr ~ scale(outmigrants.pc) + scale(inmigrants.pc)
                 + log(electorate) + log(pop.tot) + log(people_between_18_and_70) + log(people_above_16) + log(area_km2)
                 + log(pop.urb) + log(pop_higher_educ) + avg.capita.income + Gini + pc.poor 
                 | microregion + year, 
                 data = wvm, vcov_cluster(~microregion), weights = ~pop.tot)

# Table A9 
etable(m5s.amc, m3s.amc, m4s.amc, m5s.mcr, m3s.mcr, m4s.mcr,
       drop = c("pop.tot", "electorate", "people_between_18_and_70", "people_above_16", "area_km2",
                "pop.urb", "pop_higher_educ", "avg.capita.income", "Gini", "pc.poor"),
       fitstat = ~ ar2 + n, 
       dict = c("turnout_rate_fr"="Turnout Rate", 
                "scale(outmigrants.pc)" = "Out-migration", "scale(inmigrants.pc)" = "In-migration",
                "scale(turnover.pc)" = "Migration Turnover", "scale(net.migration.pc)" = "Net Migration"),
       title = "Panel Data Analysis: Turnout Rate and Migration Shares (Stock Measure), MCA and Micro-region Levels",
       style.tex = style.tex("aer"), digits = "r3", tex = TRUE)


# --------------
#   Table A10 
# --------------

# ----------- Migration Flow ----------- 

m3ow.r <- feols(turnout_rate_fr ~ scale(turnover.flow.pc) 
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                + local.election
                | id_municipio, 
                data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m4ow.r <- feols(turnout_rate_fr ~ scale(net.mig.flow.pc) 
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                + local.election
                | id_municipio, 
                data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m5ow.r <- feols(turnout_rate_fr ~ scale(outmig.flow.pc) + scale(inmig.flow.pc)
                + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                + local.election
                | id_municipio, 
                data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

# ----------- Migration Stock ----------- 

m3sow <- feols(turnout_rate_fr ~ scale(turnover.pc) 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio, 
               data = wv, vcov_cluster(~id_municipio), weights = ~pop)

m4sow <- feols(turnout_rate_fr ~ scale(net.migration.pc) 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio, 
               data = wv, vcov_cluster(~id_municipio), weights = ~pop)

m5sow <- feols(turnout_rate_fr ~ scale(outmigrants.pc) + scale(inmigrants.pc)
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | id_municipio, 
               data = wv, vcov_cluster(~id_municipio), weights = ~pop)


#  Table A10
etable(m5ow.r, m3ow.r, m4ow.r, m5sow, m3sow, m4sow,
       drop = c("pop", "electorate", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor", "local.election"),
       vcov = ~id_municipio, fitstat = ~ ar2 + n, 
       dict = c("turnout_rate_fr"="Turnout Rate", "scale(outmig.flow.pc)" = "Out-migration", 
                "scale(inmig.flow.pc)"="In-migration", "scale(turnover.flow.pc)"="Migration Turnover", 
                "scale(net.mig.flow.pc)"="Net Migration",
                "scale(outmigrants.pc)" = "Out-migration", "scale(inmigrants.pc)" = "In-migration",
                "scale(turnover.pc)" = "Migration Turnover", "scale(net.migration.pc)" = "Net Migration"),
       title = "Municipality Fixed Effects: Turnout Rate and Migration Shares (Flow and Stock Measures)",
       style.tex = style.tex("aer"), digits = "r3", tex = TRUE)


# --------------
#   Table A11 
# --------------

# ----------- Migration Flow ----------- 

m3yfe.r <- feols(turnout_rate_fr ~ scale(turnover.flow.pc) 
                 + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                 + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                 + local.election
                 | year, 
                 data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m4yfe.r <- feols(turnout_rate_fr ~ scale(net.mig.flow.pc) 
                 + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                 + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                 + local.election
                 | year, 
                 data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

m5yfe.r <- feols(turnout_rate_fr ~ scale(outmig.flow.pc) + scale(inmig.flow.pc) 
                 + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
                 + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
                 + local.election
                 | year, 
                 data = mgp, vcov_cluster(~id_municipio), weights = ~pop)

# ----------- Migration Stock ----------- 

m3yfe <- feols(turnout_rate_fr ~ scale(turnover.pc) 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | year, 
               data = wv, vcov_cluster(~id_municipio), weights = ~pop)

m4yfe <- feols(turnout_rate_fr ~ scale(net.migration.pc) 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | year, 
               data = wv, vcov_cluster(~id_municipio), weights = ~pop)

m5yfe <- feols(turnout_rate_fr ~ scale(outmigrants.pc) + scale(inmigrants.pc)
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(area_km2)
               + log(pop_urb) + log(pop.college) + avg.capita.income + Gini + pc.poor 
               | year, 
               data = wv, vcov_cluster(~id_municipio), weights = ~pop)

#  Table A11
etable(m5yfe.r, m3yfe.r, m4yfe.r, m5yfe, m3yfe, m4yfe,
       drop = c("pop", "electorate", "num_people_between_18_and_70", "num_people_above_16", "area_km2",
                "pop_urb", "pop.college", "avg.capita.income", "Gini", "pc.poor", "local.election"),
       vcov = ~id_municipio, fitstat = ~ ar2 + n, 
       dict = c("turnout_rate_fr"="Turnout Rate", "scale(outmig.flow.pc)" = "Out-migration", 
                "scale(inmig.flow.pc)"="In-migration", "scale(turnover.flow.pc)"="Migration Turnover", 
                "scale(net.mig.flow.pc)"="Net Migration",
                "scale(outmigrants.pc)" = "Out-migration", "scale(inmigrants.pc)" = "In-migration",
                "scale(turnover.pc)" = "Migration Turnover", "scale(net.migration.pc)" = "Net Migration"),
       title = "Year Fixed Effects: Turnout Rate and Migration Shares (Flow and Stock Measures)",
       style.tex = style.tex("aer"), digits = "r3", tex = TRUE)


# --------------
#   Table A12
# --------------

# Transform turnout rate
mgp$turnout_rate_fr.frac <- mgp$turnout_rate_fr/100

# Declare panel data for WBRE
panel_wbre <- panel_data(mgp, id = id_municipio, wave = year) 

wb1 <- wbgee(turnout_rate_fr.frac ~ outmig.flow.pc 
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) 
             + log(num_people_above_16) + log(area_km2) + log(pop_urb) 
             + log(pop.college+1) + avg.capita.income + Gini + pc.poor
             , data = panel_wbre, model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1") 

wb2 <- wbgee(turnout_rate_fr.frac ~ inmig.flow.pc 
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) 
             + log(num_people_above_16) + log(area_km2) + log(pop_urb) 
             + log(pop.college+1) + avg.capita.income + Gini + pc.poor
             , data = panel_wbre, model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1") 

wb3 <- wbgee(turnout_rate_fr.frac ~ turnover.flow.pc 
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) 
             + log(num_people_above_16) + log(area_km2) + log(pop_urb) 
             + log(pop.college+1) + avg.capita.income + Gini + pc.poor
             , data = panel_wbre, model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1") 

wb4 <- wbgee(turnout_rate_fr.frac ~ net.mig.flow.pc 
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) 
             + log(num_people_above_16) + log(area_km2) + log(pop_urb) 
             + log(pop.college+1) + avg.capita.income + Gini + pc.poor
             , data = panel_wbre, model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1") 

wb5 <- wbgee(turnout_rate_fr.frac ~ outmig.flow.pc + inmig.flow.pc
             + log(electorate) + log(pop) + log(num_people_between_18_and_70) 
             + log(num_people_above_16) + log(area_km2) + log(pop_urb) 
             + log(pop.college+1) + avg.capita.income + Gini + pc.poor
             , data = panel_wbre, model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1") 

# List models
wbre_models_main <- list(wb1, wb2, wb5, wb3, wb4)

# Labels
coeff_labels_wbgee <- 
  c("outmig.flow.pc" = "Out-migrant Share (within)",
    "inmig.flow.pc" = "In-migrant Share (within)", 
    "turnover.flow.pc" = "Migration Turnover (within)", 
    "net.mig.flow.pc" = "Net Migration (within)", 
    "imean(outmig.flow.pc)" =  "Out-migrant Share (between)",
    "imean(inmig.flow.pc)" = "In-migrant Share (between)",
    "imean(turnover.flow.pc)" =  "Migration Turnover (between)",
    "imean(net.mig.flow.pc)" =  "Net Migration (between)")

# Within-between random-effects models
options(modelsummary_format_numeric_latex = "plain")
modelsummary(wbre_models_main, 
             stars = TRUE, 
             coef_map = coeff_labels_wbgee,
             gof_omit = 'DF|Deviance|AIC|BIC|Log.Lik.',
             coef_omit = "year|id_municipio|(Intercept)|log(electorate)|log(pop)|log(num_people_between_18_and_70)|log(num_people_above_16)|log(area_km2)|log(pop_urb)|log(pop.college+1)|avg.capita.income|Gini|pc.poor",
             output = "latex")

# -------------
#   Table A13
# -------------

# Social Cohesion Data
soc <- read.csv("social-cohesion.csv")

# Merge with average census data 2000-2010
soc <- left_join(soc, datamean, by = "id_municipio")

# Define social cohesion measure
soc$socialcohesion.sd <- rowMeans(soc[c("ln.pc.CSOn.sd", "homiciderate.sd", 
                                        "TRUSTINSmrp.sd", "TRUSTSOCmrp.sd", "NATPRIDEmrp.sd")], na.rm = T)

m1t0222 <- lm(socialcohesion.sd ~ outmigrants.stock.pc 
              + state +  log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + Gini + pc.poor + pop_urb + pop.college +
              + log(area_km2) + distance_to_coast + distance_to_capital + longitude + latitude
              , data = soc)

m2t0222 <- lm(socialcohesion.sd ~ inmigrants.stock.pc 
              + state +  log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + Gini + pc.poor + pop_urb + pop.college +
              + log(area_km2) + distance_to_coast + distance_to_capital + longitude + latitude, data = soc)

m3t0222 <- lm(socialcohesion.sd ~ turnover.stock.pc 
              + state + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + Gini + pc.poor + pop_urb + pop.college +
              + log(area_km2) + distance_to_coast + distance_to_capital + longitude + latitude, data = soc)

m4t0222 <- lm(socialcohesion.sd ~ net.mig.stock.pc 
              + state + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + Gini + pc.poor + pop_urb + pop.college +
              + log(area_km2) + distance_to_coast + distance_to_capital + longitude + latitude, data = soc)

m5t0222 <- lm(socialcohesion.sd ~ outmigrants.stock.pc + inmigrants.stock.pc 
              + state + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + Gini + pc.poor + pop_urb + pop.college +
              + log(area_km2) + distance_to_coast + distance_to_capital + longitude + latitude, data = soc)

# Table A13
stargazer(m1t0222, m2t0222, m5t0222, m3t0222, m4t0222, 
          dep.var.caption = "",
          title="Migration Shares and Social Cohesion Across Brazilian Municipalities", 
          dep.var.labels=c("Social Cohesion"),
          covariate.labels=c("Out-migration", "In-migration", "Migration Turnover", "Net Migration"),
          no.space=TRUE, column.sep.width = "10pt", 
          omit = c("Constant", "pop", "electorate", "num_people_between_18_and_70", "num_people_above_16", "Gini", "pop_urb", "pop.college", "pc.poor",
                   "area_km2", "distance_to_capital", "distance_to_coast", "longitude", "latitude", "state"),
          omit.stat=c("rsq", "ser", "f"), notes = "", notes.append = F,  notes.label = "",
          star.char = c("+", "*", "**", "***"), star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          add.lines = list(c("State FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Geographic controls", "Yes", "Yes", "Yes", "Yes", "Yes")))


#----------------
#   Table A14
#----------------

# Survey data
dataIndividual <- read.csv("LAPOP-turnout.csv")

# Merge municipal-level data
dataIndividual <- left_join(dataIndividual, datamean, by = "id_municipio")

# Create migration rates
dataIndividual$outmigrants.stock.pc <- dataIndividual$outmigrants.stock / dataIndividual$pop
dataIndividual$inmigrants.stock.pc <- dataIndividual$inmigrants.stock / dataIndividual$pop
dataIndividual$turnover.stock.pc <- dataIndividual$turnover.stock / dataIndividual$pop
dataIndividual$net.mig.stock.pc <- dataIndividual$net.mig.stock / dataIndividual$pop

# Models
surv0 <- feols(voted ~ turnover.stock.pc 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) +
               + log(pop.college) + log(pop_urb) + Gini + pc.poor
               | state + year, 
               data = dataIndividual, vcov_cluster(~id_municipio))

surv1 <- feols(voted ~ addNA(nonmover) 
               + factor(female) + factor(agegroup) + factor(race) + factor(education) 
               | state + year, 
               data = dataIndividual, vcov_cluster(~id_municipio))

surv2 <- feols(voted ~ turnover.stock.pc + addNA(nonmover) 
               + factor(female) + factor(agegroup) + factor(race) + factor(education) 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(pop.college) + log(pop_urb) + Gini + pc.poor
               | state + year, 
               data = dataIndividual, vcov_cluster(~id_municipio))

surv3 <- feols(voted ~ turnover.stock.pc + addNA(nonmover) 
               + factor(female) + factor(agegroup) + factor(race) + factor(education) 
               + log(electorate) + log(pop) + log(num_people_between_18_and_70) + log(num_people_above_16) + log(pop.college) + log(pop_urb) + Gini + pc.poor
               + log(area_km2) + distance_to_capital + distance_to_coast + longitude + latitude
               | state + year, 
               data = dataIndividual, vcov_cluster(~id_municipio))

# Table A14
etable(surv1, surv0, surv2, surv3,
       dict = c("voted"="Voted in Last Election", "turnover.stock.pc" = "Municipal Turnover Share",
                "addNA(nonmover)1" = "Non-Migrant Status", "addNA(nonmover)NA" = "NA"),
       drop = c("NA", "female", "agegroup", "race", "education", "pop",
                "electorate", "num_people_between_18_and_70", "num_people_above_16", "pop.college", "pop_urb", "Gini", "pc.poor",
                "area_km2", "distance_to_capital", "distance_to_coast", "longitude", "latitude"),
       fitstat = ~ ar2 + n, 
       title = "Individual Voting Behavior, Migration Status, and Local Turnover",
       style.tex = style.tex("aer"), digits = "r3", tex = TRUE)






